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It is now possible to compute linear in mass-ratio terms in the post-Newtonian (PN) expansion for 
compact binaries to very high orders using linear black hole perturbation theory applied to various 
invariants. For instance, a computation of the redshift invariant of a point particle in a circular 
orbit about a black hole in linear perturbation theory gives the linear-in-mass-ratio portion of the 
binding energy of a circular binary with arbitrary mass ratio. This binding energy, in turn, encodes 
the system’s conservative dynamics. We give a method for extracting the analytic forms of these 
post-Newtonian coefficients from high-accuracy numerical data using experimental mathematics 
techniques, notably an integer relation algorithm. Such methods should be particularly important 
when the calculations progress to the considerably more difficult case of perturbations of the Kerr 
metric. As an example, we apply this method to the redshift invariant in Schwarzschild. Here we 
obtain analytic coefficients to 12.5PN, and higher-order terms in mixed analytic-numerical form 
to 21.5PN, including analytic forms for the complete 13.5PN coefficient, and all the logarithmic 
terms at 13PN. We have computed the individual modes to over 5000 digits, of which we use at 
most 1240 in the present calculation. At these high orders, an individual coefficient can have over 
30 terms, including a wide variety of transcendental numbers, when written out in full. We are 
still able to obtain analytic forms for such coefficients from the numerical data through a careful 
study of the structure of the expansion. The structure we find also allows us to predict certain 
“leading logarithm”-type contributions to all orders. The additional terms in the expansion we 
obtain improve the accuracy of the PN series for the redshift observable, even in the very strong- 
held regime inside the innermost stable circular orbit, particularly when combined with exponential 
resummation. 

PACS numbers: 04.25.Nx, 04.30.Db, 04.70.Bw 

I. INTRODUCTION AND SUMMARY 

Coalescing compact binaries are a promising source 
of gravitational waves, and ground-based gravitational 
wave interferometers will start operating at sensitivities 
at which detections can reasonably be expected as early 
as later this year. In order to successfully detect these 
faint signals in the detector’s noise—and, more impor¬ 
tantly, to be able to infer the properties of the system 
from the detected signal—it is necessary to have highly 
accurate templates that model the gravitational waves 
from the inspirailing binaries. Thus, for more than a 
decade now, different approaches have been developed 
to model relativistic binary systems. The oldest one of 
these, the post-Newtonian (PN) framework, can model 
such systems when the two bodies are far from one an¬ 
other, so their velocities are relatively slow (see [I| for 
a review of these methods and results). Numerical rel¬ 
ativity, on the other hand, is able to model compara¬ 
ble mass ratio binaries in the strong gravitational held 
regime, but has difficulties with large mass ratios, large 
separations, and very long waveforms (but see 
recent advances). Another approach is gravitational self¬ 
force theory, which models binaries with extreme mass- 
ratios, where one has a small body that is about a million 
times lighter than the central super-massive black hole 


into which it is inspiralling [H3- 

A more recent approach, effective-one-body (FOB) 
theory, maps the binary’s motion to that of a particle 
moving in an effective metric, generalizing the Newtonian 
reduced-mass treatment of the two-body problem [1, . 

This theory encompasses information from the former 
three approaches to calibrate the parameters that go into 
the theory, which allows it to model a binary system of 
any given mass-ratio. Of particular interest is the over¬ 
lap region between the self-force and PN formalisms. In¬ 
variant quantities calculated in this region are used to 
calibrate the EOB parameters. One of those quantities 
calculated in self-force theory is Detweiler’s redshift in¬ 
variant, A[/, the linear-in-mass-ratio correction to the 
time component of the 4-velocity of the light compact 
object [^. The PN coefficients of AC7 are directly re¬ 
lated to those of the linear-in-mass-ratio portion of the 
binding energy and angular momentum of the binary, as 
well as to the radial potential that is fundamental to the 
EOB formalism, as was demonstrated in [IMl- 

Computation of the PN coefficients of AN started with 
Detweiler’s original paper (to 2PN; nPN corresponds 
to an accuracy of where v is the orbital velocity of 
the small body), and continued analytically through 3PN 
in 0: using standard post-Newtonian methods, with 
terms through 5PN obtained from a numerical match- 
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iiig in [l^ (where the logarithmic terms were obtained 
analytically). In [l^, Bini and Damour calculated the 
nonlogarithmic portion 4PN coefficient analytically by 
using analytical solutions of the Regge-Wheeler-Zerilli 
equations; the logarithmic term had already been com¬ 
puted by Damour in [ 2 ^. Shah, Friedman, and Whit¬ 
ing [HI (hereafter referred to as SFW) calculated higher 
order PN coefficients, up to 10.5PN order, by calculat¬ 
ing high-precision numerical values in a modified radi¬ 
ation gauge at very large radii and fitting them to a 
PN series to extract the coefficients. They found that 
half-integer terms started at 5.5PN, which they veri¬ 
fied analytically (and was also verified using standard 
PN methods in j^, [ 2 ^). SFW were also able to in¬ 
fer analytic forms for certain not-too-complicated coeffi¬ 
cients from their numerical data. Concurrently, Bini and 
Damour calculated the coefficients analytically to 6PN 
order in M- Analytical calculations since then [HI - It 7| . 
all using the Regge-Wheeler-Zerilli gauge and the results 
from |28| . have been pushed to find much higher order 
PN coefficients. We have compared our results to 20.5PN 
with the concurrent calculation by Kavanagh, Ottewill, 
and Warden , who have streamlined the Bini-Damour 
method, and found complete agreement. (Note that Bini 
and Damour’s 9.5PN result appeared while we were 
finishing checking higher-order terms in this work.) 

Apart from AU, high-order PN coefficients of other in¬ 
variants from which the FOB formalism can benefit have 
also been calculated: These are the linear-in-mass ratio 
conservative corrections to the spin-precession angle [ 2 ^ 
and the quadrupolar [sO] and octupolar [H| tidal invari¬ 
ants (the eigenvalues of the electric- and magnetic-type 
tidal tensors)jall which have been calculated to high 
PN order in jH> Hfj-fU. Recently, Bernuzzi et al. 
introduced a semi-analytical tidally coupled binary neu¬ 
tron star model using the FOB theory where information 
from the PN expansion of the redshift and tidal invari¬ 
ants was incorporated (using the tidal results from Bini 
and Damour [^). The results of this model are in good 
agreement with a more recent numerical simulation in full 
general relativity by the Japanese school (Hotokezaka et 
al. [ 3 ^) in the case of compact neutron stars. 

It has recently been shown (see maiig) how the 
overlap region between the self-force formalism and the 
PN approximation can be explored using very high ac¬ 
curacy numerical results, which make it relatively easy 
to extract high-order PN coefficients that are currently 
out of reach of standard PN calculations. The coeffi¬ 
cients obtained using this numerical extraction method 
have then been checked by independent analytical calcu¬ 
lations. The advantage of developing such high-accuracy 
calculations will be evident when PN coefficients are cal¬ 
culated for invariants in Kerr spacetime, where purely 
analytical calculations will likely be extremely difficult. 
The techniques developed in this paper can then be gen¬ 
eralized to calculate the analytical form of the numerical 
coefficients for various invariants in Kerr, and, eventu¬ 
ally, to calculate the quadratic-in-mass-ratio terms us¬ 


ing second-order self-force results (see Sec. 4.3.3 in Q 
for a brief overview of progress in second-order self-force 
calculations). Additionally, one obtains insight into the 
structure of the PN expansion from these high-order com¬ 
putations, particularly in comparing the terms one can 
predict using a simplification in AC/ with similar terms 
in the energy flux at infinity [38| . 

We shall now outline our method and compare it to 
previous work. Shah, Friedman, and Whiting (SFW) [2l| 
worked solely on the expansion of the full AU and ob¬ 
tained analytic terms for the simplest coefficients, which 
are purely rational, or a rational times tt, where they 
could easily identify the analytic form from a large 
enough number of digits. They also present three addi¬ 
tional analytic expressions for more complicated higher- 
order terms (in the note added), which were obtained by 
the first author of this paper using an integer relation 
algorithm. However, the accuracy of the expressions in 
SFW was insufficient to obtain analytic forms for any 
terms beyond 10.5PN order. 


The methods we use here are similar to those used to 
obtain the more complicated coefficien ts g iven in SFW 
(and the analytic coefficients given in 122 1), in that 

we also use an integer relation algorithm, but the present 
application is more effective at obtaining higher-order 
terms, since we primarily work with the individual modes 
of AU [either retarded {£, m) modes, or renormalized £ 
modes], where the structure of the expansion is simpler, 
and one can obtain analytic forms at a given order with 
fewer digits. Indeed, one can often predict some—and 
in certain cases even all—of the entire analytic form at 
higher orders from lower-order coefficients. This simpli¬ 
fication of the structure when considering the individual 
(£, m) modes was also seen in the expansion of the en¬ 
ergy flux at infinity of a point particle in a circular orbit 
around a Schwarzschild black hole [1^ . Additionally, the 
overall structure of the expansion of the retarded (£, m) 
modes of AU is also similar to that of the energy flux 
at infinity (calculated to 22PN by Fujita [s^, with the 
structure studied in [l^), and we are able to use this to 
help determine which transcendentals to include in the 
vector to which we applied the integer relation algorithm. 

We also use the integer relation algorithm in a more 
fundamental way in the current work, preferring for most 
of our work to find analytic expressions for the terms 
order-by-order and then subtract them off to obtain the 
numerical values for higher-order terms to higher accu¬ 
racy. (Note that we found that for certain of the more 
complicated terms at higher orders it was necessary to 
first obtain analytic forms for some of the simpler co¬ 
efficients at even higher orders in order to obtain the 
more complicated terms to sufficient accuracy to be able 
to determine an analytic form.) This method should be 
contrasted with the more usual method of finding nu¬ 
merical values for all terms to some accuracy using a fit, 
then finding analytic forms for some coefficients, using 
these to improve the accuracy of the fit, and iterating. 
This fitting method was used in SFW and in Nickel’s 
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similar computation of high-order terms in the expan¬ 
sion of t he g round state energy of H 2 in powers of the 
distance [4^.^ We also used this fitting method on the 
full AU to obtain even higher-order terms than we were 
able to obtain using the first method, though these terms 
were all obtained only in mixed numerical-analytic form. 

These i nteg er relation algorithms, notably the PSLQ 
algorithm |42l . l43l | , are a prominent tool in modern exper¬ 
imental mathematics. (See also 0 1 ^ for further intu¬ 
ition into the PSLQ algorithm and for a review of 
some remarkable results obtained using integer relation 
algorithms. Additionally, see [Pzl - ld^ for some general 
reviews of the methods, philosophy, and results of mod¬ 
ern experimental mathematics.) The PSLQ algorithm 
returns a small vector of integers that is orthogonal to a 
given input vector, and thus can be used to identify the 
analytic form of numbers from a high-accuracy decimal 
expansion, which is the task for which we use it here, em¬ 
ploying the implementation in the FindIntegerNullVector 
function in Mathematica (first available in version 8). 

Here we only need to identify numbers that are lin¬ 
ear combinations of transcendentals with rational coef¬ 
ficients, which is one of the simplest cases to which one 
can imagine applying an integer relation algorithm. Nev¬ 
ertheless, there are enough transcendentals at higher or¬ 
ders, with complicated enough rational coefficients, that 
we still need to compute certain individual PN coeffi¬ 
cients to over 200 digits, even when using a simplification 
we found that helps remove much of the complexity at 
higher orders. This necessitates computing the modes of 
AU to over 1000 digits; we actually computed to over 
5000 digits so we could go to even higher orders, where 
we currently only obtain certain coefficients analytically. 
Some other such high-accuracy computations in mathe¬ 
matical physics, including further applications of PSLQ, 
are discussed in Additionally, Nickel 0 also uses 
PSLQ to obtain analytic forms of high-order coefficients 
of a similar series for the ground-state energy of 

We also note, following Bini and Damour [^, that 
while one has to sum over all spherical harmonic (£, m) 
modes to obtain AU to a given PN order (compared to, 
e.g., the energy flux, where one only has to sum a finite 
number of modes to obtain the expansion to a given PN 
order), this infinite sum is only necessary to obtain the 
nonlogarithmic integer-order PN terms. All the other 
terms in the PN expansion of AC/ come from a finite 
sum over modes. 

It also turns out that the expression for the PN coef¬ 
ficient of a given renormalized i mode (at high £, where 
it is purely rational) is simple enough that one can infer 
it from the numerical values of fewer than 100 C-modes, 


^ This expansion of the ground state energy of has a similar 
structure to an individual mode of AU (though it is simpler), 
and can be computed using functional series techniques similar 
to those we employ here, as discussed in |41| . 


at the PN orders at which we are working, so we can 
obtain these general expressions and then perform the 
infinite £-sum analytically, allowing us to calculate ana¬ 
lytic forms for the nonlogarithmic integer-order PN coef¬ 
ficients of AU without performing the Asum numerically. 
This is fortunate, since performing the infinite /-sum nu¬ 
merically to such high accuracies would be prohibitively 
expensive, computationally, due to the necessity of calcu¬ 
lating many i modes. We obtained the full expansion up 
to 12.5PN this way (including reproducing all the known 
analytic terms “from scratch”). 

We can even obtain fairly complicated forms of high- 
order terms (though not any complete PN terms) using 
the predictions of the simplification, PSLQ, and a rea¬ 
sonably (but not excessively) accurate calculation of the 
full AU. In particular, we performed a calculation of the 
full AU to “merely” ^ 600 digits at somewhat smaller 
radii (10^®M to 9 x lO^^M, where M is the mass of the 
central object) to check the values we obtained using the 
data calculated to more than 5000 digits for radii from 
10®°M to 10™M (where we used at most 1240 digits and 
15 radii to obtain those results).^ Using the results of 
this calculation, we were able to obtain accurate values 
up to 21.5PN, including analytic forms for 48 coefficients 
containing as many as 27 terms (most coefficients had 
far fewer terms), starting from the full AU, though most 
of these terms were predicted by the simplification: We 
only had to use PSLQ to obtain at most 4 terms. Here 
we also used the analytic forms we obtained to iteratively 
improve the PN coefficients, increasing the accuracy of 
terms we had already obtained, in addition to obtaining 
even higher-order terms. 

The plan of the paper is as follows. We first briefly 
review the relevant portions of the self-force calculation 
in Sec. m and then discuss the method we use to ob¬ 
tain the PN coefficients of the individual modes of AU 
(including a simplification of the modes, and consistency 
checks) in Sec. lHIl We then give the terms in the full AU 
that are predicted to all orders by the simplification of 
the modes in Sec. HYi and discuss how we compute the 
infinite sum over the modes of AU to obtain the final re¬ 
sults for the PN coefficients of AU (and our independent 
check of these results) in Sec. |Vl We discuss convergence 
of the series in Sec. ivn and conclude in Sec. IVHI In the 
Appendix, we give some discussion of how one can obtain 
certain parts of the simplifications of the modes of AU 
from inspection of the method we use to calculate it. We 
use geometrized units throughout (setting the speed of 
light and Newton’s gravitational constant both to unity, 
i.e., G = c = 1). 


^ These calculations were computationally not exceptionally ex¬ 
pensive, requiring ~ 45 hours per radius on 2 processors at ~ 600 
digits, and ~ 10 hours per radius on 16 processors at ~ 5000 dig¬ 
its. 
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II. SELF-FORCE CALCULATION 


components, u* and are given by 


Here we give the basics of the method we use to calcu¬ 
late At/ (and its precise definition)—see [U for 

further details. We calculate AC/ in a modified radiation 
gauge, where / > 2 modes are calculated in an outgoing 
radiation gauge (with hapn'^ = 0 and h = 0, where n“ 
is the ingoing null vector and hap and h are the metric 
perturbation and its trace, respectively) and the lower 
ones (/ = 0,1) are calculated in the asymptotically flat 
Schwarzschild gauge. The setup is as follows: A particle 
of mass m is orbiting a Schwarzschild black hole of mass 
M in a circular orbit of radius r = tq in Schwarzschild 
coordinates (/, r, 0, cj)). The particle’s four-velocity, u“, is 
given by 


=uH°‘ (1) 

where and <()“ are the time-like and rotational Killing 
vectors of the Schwarzschild metric, respectively. The 


U :=u* = 



with 


n = 



( 2 ) 

( 3 ) 

( 4 ) 


We follow the Chrzanowski-Cohen-Kegeles-Wald formal¬ 
ism (outlined in [5l|) of extracting the metric pertur¬ 
bation from the perturbed spin-2 Weyl scalar i/jq as fol¬ 
lows. We first solve the spin-2 separable Teukolsky equa¬ 
tion, whose retarded solution, tpQ (the superscript “ret” 
is omitted here), is given by 


i’oix) = ' 00 °^ + ^ 0 ^^ + ( 5 ) 


with 


^0 im ^ 

= 87riTnHu‘Ao ^ - 1)(^ + 2)]^/^2Wm(6', (j))lYtra ^ 

Em 

|[iTOr2ro -I- 2ro]i?H(c<)i?oo(r'>) -I- Ao[R^{ro)R^{r)e{r - tq) -|- R}i{r)R'^{ro)0{ro - r)]|, (6b) 

= -4:Trmfl^u*'^Aem2Yem(0,<P)2yem(^-^,^t^x 

Em 

|[30rQ — SOMtq + LSM^Cq — — 2Aq — 24Aoro(ro — M) + 6im^lrQ{rQ — M)]Rii{r^)Roo{ry) 

-f2(6r^ - 20Mrp -H - SrpAg -h imHAorp)[i?H(ro)i?oo(c)6»(r - tq) -/ R'^{ro)Riiir)0iro - r)] 

+rQAl[R'^{ro)Rao{r)0{r - tq) + R'^{ro)R}iir)0{ro - r) + W[i?H(r),i?oo(?')]^(?' - 'Co)] i, (6c) 


where A := — 2Mr; the function Rh is the solution 

of the homogenous radial Teukolsky equation which is 
ingoing at the future event horizon and i?oo is the one 
that is outgoing at null infinity. Here r< := min(r, tq) 
and r> := max(r, rg). A prime denotes a derivative 
with respect to the r-coordinate and overbars denote 
complex conjugation; S and 0 denote the Dirac delta 
distribution and Heaviside theta function, respectively. 
The Wronskian of these two retarded radial solutions is 


W[AH(r),i?oo(r)] = RnR'ac - The quantity 

given by 

A3W[i?H(r),i?oo(r)]’ 

is a constant independent of r, that is, A'f^^ = 0. The 
functions Rr and i?oo are calculated to more than 5000 
digits of accuracy using the Mano, Suzuki, and Takasugi 
(MST) method given in [s^, namely 
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OO 

Rh = ^ anF{n + + I — i€,—n — — ie,—l — 2ie;x), 

n— — oo 

OO 

Roo = F^z''~^ ^ (—2z)”6„J7(n + + 3 — ie, 2n + 2^ + 2; —2iz). 

n— — oo 


(8a) 

(8b) 


where a; = 1 — e = 2Mm^l, z = —ex, F is the hy¬ 
pergeometric function 2 F 1 , and U is the (Tricomi) con¬ 
fluent hypergeometric function. To expedite the calcu¬ 
lation, and reach the high accuracies we require, we use 
various recurrence relations for U, and Gauss’s relations 
for contiguous functions for 2 F 1 (see, e.g., M) to write 
the various n-dependent functions and their derivatives 
in terms of the functions calculated for n = 0,1. For de¬ 
tails regarding the derivation of v, the renormalized an¬ 
gular momentum, and the coefficients a„ and please 
refer to [13, The spin-weighted spherical harmon¬ 
ics sYiTn{0, 4>) s-re calculated analytically and are given in 

From 'ijjQ, we compute the intermediate Hertz potential, 
'h, from which the components of the metric perturbation 
are calculated. The radial parts of and '0o are related 
by an algebraic relation given by 

O (-1)'"(^ + 2)(^ + IW - l)^^,-m + 12imMHV-^™ 

[{£ -f 2)(£ -f !)£{£ - 1)]2 -p 

(9) 

where 

^ = (10a) 

i,m 

^o = Yl 0)e-’’"^*. (10b) 

l,m 


Once we compute dr, the components of the metric per¬ 
turbation along the Kinnersley tetrad are given by 


hll = y(3"^ + S"d'), 


/i.s.s = r'^ 


+ 


dj - 2fdtdr + pdl 3(r - M) 


4 2r2 

/(3r-2M)^ , 7-2-2M2 


dt 


2r2 




( 11 ) 




The linear-in-mass-ratio correction to the time compo¬ 
nent of the four-velocity of the particle due to its finite 
but small mass is then given by 


AG = (13) 


where 


^ren 




(14) 


The super-script “ren” denotes the renormalized, 
singularity-free part of the metric perturbation. We refer 
the reader to [U, Isil - f^ IstI - Is^ for details pertaining to 
the renormalization procedure. 

As mentioned earlier, ^0 only provides us with the ra¬ 
diative part of the metric perturbation. One also has to 
add on the non-radiative parts associated with the change 
in mass and angular momentum of the Schwarzschild 
spacetime with the particle. These contributions are 
given by (Eqs. (137) and (138) of [Hi) 


Hsm 

Hsj 


m(ro — 2M) 
rl^‘^{ro — 3M)3/2 
—2Mm 

rl^^{ro — 3M)3/2 


(15a) 

(15b) 


The index SM and SJ refer to the parts coming from 
the change in mass and angular momentum, respectively. 
Also note that the £ = l,m = ±1 (even) contribution 
to such gauge-invariant quantities, corresponding to the 
shift in center-of-mass of the binary m — M system, is 
zero. 


III. OBTAINING ANALYTIC FORMS OF THE 
PN COEFFICIENTS OF THE INDIVIDUAL (£, m) 
MODES OF AG 

A. The PN expansion of the (2, 2) mode of AG and 
a way to simplify a general (I, m) mode of AG 


where / = A/r^, and the angular operators S and 3, the 
s-raising and -lowering operators, are given by 

Br] = — {dg + i esc 9d^ — s cot 9) rj, 

Bt] = — {dg — icsc9d^ + scot9)ri, ^ ^ 


We start by giving our expression for the PN expan¬ 
sion of T 22 , the (2,2) mode of AU/U, through 12.5PN, 
as well as the simplification of the modes we have discov¬ 
ered, before describing our method for obtaining these 
results. [We give the analogous results for the other 
modes in the electronic Supplemental Material (H : along 
with higher-order PN coefficients in the (2, 2) mode for 


where ry has spin-weight s. 
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which we only know analytic forms for some of the terms, 
and the 13 . 5 PN piece we do know all of.] Here we con¬ 
sider AUjU instead of just At/ as this is the quan¬ 
tity that we worked with on the level of the individual 


modes [note that AUjU = —cf. Eq. (fT^ j. We 
present the expansion in terms of the same dimension¬ 
less and gauge invariant radius variable used in SEW, 
viz., R := = tq/M'. 


T 22 = 


31 1 I91 1 


1 


5813077 ^ 1 


22 R 2371772 253271 773 26315272111774 

29 I 71 I 733 I 33 OII 33231 261511 


1917271789381591 2^ , , 

29335273111131 - ^ eulerlog 2 ( 77 ) 
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-b 
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335473 


eulerlog2(i?) 


7712.5 


-b 


-b 


43 I 193 I 373 O 88739 I 6 II 


3*5375131 

I493II855571 


325372 '.V”/ ■ 365574111 335372 

30942707693267131529565293293479262815467669462411 

2243115777114133173192232291 

2323 I 2459 I 5245709091 , , 21268991 , , 

—C(3)- 03.1^2 C(5) 


TT^ ) eulerlog2(-R) 


log(2/R) 


21315372 1 3152 

I 9 I 37 I 227 I 4 OII 92033 I 72758215409831 
2*395476112132171 


29 l 07 i 6899 i_^ 

^ 355373 365374131 

2I6 2*31107! 

log(2/i?) eulerlog2(R) -b — log^(2/R) 


335172 


29107208991 

365274 


5271 

212107 I 68991 


7 r 2 - 


345273 


C( 3 ) j eulerlog2(i?) 
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224723OO973M8I99671 


375475131 

2^802849557341891 


eulerlog2 (i?) + 

214107 I 1511 


2 IIIO 7268991 

365374 


eulerlog2 (R) 


385373111131 

21410721511 




345372 


345373 


7 r 2 j eulerlog2 (-R) + 


24310731511 

345474 


335373 


eulerlog2 (i?) 


1 , [5623204122099819885533688671 

^ [ 243105576112132171191 

24107I1213I19I63I1357042731 


7713.5 


O 


385575111131 

(i^)- 


(16) 


Here 

eulerlog^(i?) := 7 + log(2m/i?i^^), (17) 

where 7 is the Euler-Mascheroni gamma constant, is the 
function associated with many higher-order tail terms in 
the PN expansion, first introduced in general by Damour, 
Iyer, and Nagar [^, with a slightly different definition, 
since they use a different variable. Additionally, ( de¬ 
notes the Riemann zeta function. 


The PN expansion of the (2, 2) mode for At/ has quite 
a bit of structure that is readily apparent in its prime fac¬ 
torization, and the PN expansions of the other modes dis¬ 
play similar structure. In particular, we can write most 
of the eulerlog()j(i?), half-integer, and zeta function terms 
(including the even powers of tt) in [the {£, m) mode 
of At//t/] to the orders currently known in the following 
form (“C” is for “complications”) 


J 


'Y'Cl _ 

^ im 


eulerlog, 


r,3 


.{R)[ 


1 1 


\ 12 i?2 




5 ^ 7^'2^C(3) ,2C(3) 

1 


• -.2 

Iji 

3 i?2 


^4C(5) 


24m RO .5 


eulerlog, 


mvem 7r3 7rC(3) , tt® 

■ - 


3 R3.5 


d«)/ 


\ 2i/ejn R'^ 




R3.5 

^C(2) - 


1 


45 i?6-5 J 29im R^ 


15 

n 00 . (fc) 

E 


— {2m)'' 




R 8 4 m RO -5 






J^k+i+l+Efr. 


1 


( 2 m)^ C(3) - ^7'^mC(4) 


1 ,4C(5) 


^ [1 - 7,»C(2)] [(;(2) + f,„C(3)1 ^ 

(fc) 


1 1 
2p£rn 7?2 


E 

A;=0 


a(^) 




00 . 

[ 1 ] 


— rW V 


fc =0 


R^ ' 


(18) 


Here we have given two forms for to better illustrate 
its structure;3 recall that ((2) = and C(4) = 7r4/90. 
Additionally, 


= (19) 


^ Note that we only need to use to simplify the m ^ 0 modes. 
The m = 0 modes are nonradiative, and thus already have purely 
rational simple integer-order PN series, with no simplification 
necessary. Therefore, even though = 0 for m = 0, one does 
not need to be concerned about potential division by zero [or 
logarithms of zero in eulerlog^(i?)] in Eq. (I18II . 


where u is the renormalized angular momentum intro- 
duced in the MST formalism [sj, ■ (Here we denote 
its dependence on £ and m explicitly, which is usually 
not done in the literature, though we suppress its de¬ 
pendence on i?, even though we displayed the analo¬ 
gous dependence on v in [38|.) See the Appendix of 
Bini and Damour for explicit expressions for 
k G {1,2,3}, where these are referred to as V 2 k{£)- Note 
also that [ 1 / 2 ], = —IO 7 I/ 2 I 3 I 5 I 71 , which explains the ap¬ 
pearance of factors of 107 in many places in the prime 
factorization of T 22 in Eq. (flHl) . In fact, v (along with 
its analogue for £ —£ — 1) gives many of the leading 

logarithms in the homogeneous solutions of the Regge- 
Wheeler equation, as noted in Sec. H B of [l^ (some 
similar results for the Teukolsky equation are also im¬ 
plicit in the results of [s^). Additionally, —iv is the 

















































monodromy of the radial Teukolsky equation about the 
irregular singular point at infinity, as is mentioned in . 
One also sees multiplied by a rational with small 
prime factors, appearing in the coefficients of integrals in¬ 
volving an i multipole in the standard PN calculation of 
the next-to-leading two half-integer terms in AU in [ 2 ^ ; 
cf. their Eqs. (3.14)-(3.18) and (4.7)-(4.11) with the val¬ 
ues for i G {2,3,4} given in Table I in [s^. Fi¬ 
nally, the general form of appears in the coefficient 
of log To (where tq is the constant associated with the 
regularization parameter tq, not to be confused with the 
To in this paper) in post-Newtonian expressions for all 
the mass-type radiative multipole moments; cf. Eq. (3.9) 
in and Eq. (A2) in [ 2 ^. This coefficient was de¬ 
rived by Blanchet and Damour in the Appendix of [g^ 
using methods that differ from both the continued frac¬ 
tion method of MST [sj and the monodromy method of 
Castro et al. [g^. We also define 

0 if ^ -I- m is even, , , 

1 if £ -I- m is odd. 


The coefficients are rational and are given by the 
coefficients of the eulerlog^(i?) terms. While we might 
expect there to be contributions to the eulerlog terms 
that are not part of this simplification starting at 9PN, 
by analogy with the remainder of the Sim factorization 
of the modes of the energy flux from [s^ , it appears that 
this is not the case, since we see the same structure in 
the remainder with this choice for the coefficients as 
for the Sim factorization of the modes of the energy flux 
to all the orders we have considered. 


If we apply this simplification to the (2, 2) mode, then 
we have 


CXD 


E 


4 W 

^22 

R>^ 


2 ^ ^ 261511 1 , 24G8991 1 2315291219111 1 , 53i4G7003G811Gli 1 , 157i3745GGlG009i 1 
^ 315171 R 335172 345372111 + 365373111131 + 21365274131 ^ 

23I37I179II871I7907I58OO78G71 1 , 179il229384882345812G3017407i 1 
23365574112131171 + 24395676112132171191 W 

I 9 I 37 I 227 I 4 OII 92033 I 72758215409831 1 

28395476112132171 + ■■■ 

( 21 ) 


which yields 


T22-T 


Cl 




E 1122 


fc=i 


-k 


26 1 


2632 1 2645071 1 22191I1243431 1 


I 493 II 855571 1 


51 i?io 315171 All 


325372 7712 


21315372 A13 


log(2/A) 


246 log(2/A)eulerlog2(A) 

■261074 1 

2*341074 1 

log2(2/A) + 

24OIO74 1 

3452 

Ai3 


345274 A42 

5274 A13 

325274 A42 

2911 I 132 1 

7 r 2 -h 

r 244 7 

24032 1 

243 

TT 242194 

TT 240111 G 734 TT 

325274 7713 

54 A12 

54 A13 

3452 A 41-5 3252717712.5 345271 7713.5 


-k 


/J- 

1^14 

2 I 6 IO 7 I 

315371 


and higher terms that we do not yet know all of 


29240013G371 2ii'l07i 


)17 


log(2/A) - ^TT 


TT 2131072 log3(2/A) 


A144 


335372 A16 


o 


335472111 
1 


325371 


A164 


eulerlog 2 (R) 


( 22 ) 


(k) 

where 0 : 22 ’^ ^ Q- Here, to simplify the remainder, we 
have not used [v]^ in ^ 22 , but rather the expression for 
the l/A® piece of the 1/A expansion of ly that is valid 
for all 7 > 2, without the additional piece that only con¬ 
tributes for £ = 2 (for positive i). Specifically, this is the 
expression for i^ei^) given in the Appendix of Bini and 


Damour with the final cq term omitted. We omit 
this term because cq = 2174/3^5^1074, and we do not ac¬ 
tually see such factors of 107 (or any other anomalously 
large primes) in the denominators of the PN expansion 
of fluxes or gauge-invariant self-force quantities available 
to date. For instance, if we had used [v]^ instead of the 
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general (. > 2 expression in T 22 , then the 11.5PN and 
higher half-iirteger coefficients in the remainder would 
have had more complicated expressions (with factors of 
107 in the denominator). In particular, the 11.5PN coef¬ 
ficient in the remainder would be 2^^6737^7r/3^5^7^107^. 
However, we do see such factors of 107, and other large 
primes from the numerator of the in the denomina¬ 
tors of certain terms in the factorizations of the modes 
of the energy flux at infinity given iir [3^ . which is not 
surprising, since the full iz is used in these simplifications. 

Indeed, if one looks at the coefficients of eulerlog^(i?) 
in the PN expansion of the logarithms of the modes of 
the energy flux at infinity (for a point particle in a circu¬ 
lar orbit around a Schwarzschild black hole), then these 
coefficients give the coefficients of the 1/R expansion of 
v up to the point at which the first departure from the 
general ^-behavior occurs (at Specifically, if 

rjim denotes the {£, m) mode of this energy flux, the co¬ 
efficient of eulerlog^(r;)u®” in logTy^m is 2 [vim]n (2m)^” 
for n < £ -\- 1. This behavior is to be expected, given 
the action of the Sim factorization from (38j| ; n = £ -\- 1 
is the point at which the —v — 1 portion of 7?^“^ in the 
MST formalism starts to contribute eulerlog^ (r;) terms 
to rjim [cf. Eq. (19b) in and the discussion in Sec. IV 
of that paper]. What is striking is that the coefficient 
of eulerlog^(r')r;®(^+^i is given by the general-^ expres¬ 
sion for the PN coefficients of v. For an illustration of 
all of this for the (2, 2) mode, compare the expression 
for log?7^m in Eq. (32) in [s^ with the expressions for 
the PN expansion of u given in the Appendix of Bini and 
Damour [24| . However, for n = £-1-2 things are more com¬ 
plicated, since the expression for the PN coefficient for a 
general £ develops a pole, whose residue seems to have 
nothing to do with the difference between the coefficients 
of eulerlog^(r;)u®” in logTy^m and the PN coefficients of 
V. In particular, this residue is quite simple and has no 
large prime factors. 

The simplihcation we introduce here is in many ways 
analogous to the S(rn factorization of the modes of the 
energy flux at infinity introduced in [ 3 ^ . and likely has 
a similar expression in terms of gamma functions, where 
the current expression is just low-order terms in its PN 
expansion. However, while it is reasonably easy to read 
off the Sem factorization from the MST expression for the 
modes of the energy flux, it is far less easy to ascertain the 
similar full expression for this simplification of the modes 
of AC/, except for the euieriog„(i?) g^g discussed 

in the Appendix. Note also that the Sem factorization is 
applied to the entire mode (by division), while here we 
only subtract off a portion of the expansion with the 
simplification. 

There is also some notable structure in the remainder 
(in particular all the factors of 107 in the prime factor¬ 
ization), and it appears that the powers of log(2/i?) and 
the (^(3) terms in the remainder can all be derived from 
a single series, akin to T^. By analogy with T 22 and 
the Vf.m. factorization of the modes of the energy flux 
from (33 (though the Vim factorization does not remove 
some of the terms that are analogous to those considered 


here, so the analogy is far from exact), we conjecture that 
it has the form 


"rC2 _ I 2iyt„ log(2/ii) 
^ Im ~ ^ 


1 


2 z^£y] 


E 


B 


(k) 

im 


^A;+ 5 + 2 ^+€^t, 


TDK 

[ 2 ] ^Im 


_. r^[ 2 ] 

— ■ '-'im 2^ 


k=0 




^3 


1 


2 r^^rj 


(23) 


where 

26 2^32 1 2345073 1 

¥ ~ ~^R ~ 315171 R ^ 

22191I1243431 1 I493II855571 1 

325372 R 3 21315372 R ^ ■ 

( 24 ) 

and we give the expressions for the other modes to the 
order we know them in the electronic Supplemental Ma¬ 
terial [o^- However, note that we do not yet know the 
expansion of the individual modes to high enough orders 
to be able to check whether many of the predicted terms 
appear, and whether the coefficient of log(2/i?)/i?i^ also 
gives the coefficients of the other higher-order terms this 
expression suggests it will. Nevertheless, we are able to 
check some of these predictions for the first appearance 
of a given power of a logarithm in the (2, 2) mode us¬ 
ing our results for the higher-order PN coefficients of the 
full AU and the rest of the simplihcation, as discussed in 
Sec. IV Al Additionally, we obtain a few less direct checks 
on more of these predictions for the (2,2) mode and oth¬ 
ers from the simplihcations of the remainders of other 
logarithmic terms in the full AU. Moreover, there is a 
very similar structure in the remainder of the Sim fac¬ 
torization of the modes of the energy hux at inhnity , 
lending further support to this conjectured form. In gen¬ 
eral, these similarities between the structures that can be 
simplihed for the energy hux and AU are not surprising, 
since they likely all come from tail effects. 


00 „(fc) 

-°22 


2 ^ 

k=0 


B. Applying PSLQ to the coefficients of the PN 
expansion of the modes of AU 

We now outline the general method we use to obtain 
the analytic forms of the coefficients of the PN expansion 
of the modes of AU. First, we note that the form of the 
PN expansion of the modes of the energy hux generally 
provides a good guide to the growth of complexity of the 
terms in the expansion of the modes of AU, in particu¬ 
lar concerning the appearance of terms that we are not 
able to remove using the simplihcation—cf. the discus¬ 
sion in [ 3 ^. Next, we note that the individual retarded 
(£, to) modes of AU have the following general structure 
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E 




( 0 ) 


n—l-\-eirr 


+ TT 


i?" 


E 


E 

n—3+^+efr 

d( 0 ) 

J^r). 


eulerlog^(i?) 


i?" 


^ A^n^ eulerlog^(j?) 

n=6+e+eim 


R” 


Jln+l/2 

.n—4:-\-i-\-€^rn n—7+-^+e^ 


eulerlog^(j^) 


Bn^ eulerlog^(i?) 


+ E 

n=5+2^+e^. 


log(2/i?) 


= E 


a: 


R 

(0) oo - 


+ E 

n — ^-\-‘2^-\-sg_r 


log'(2/i?) 


R" 


n—10-\-£-\-eer} 


+ 


Jln+l/2 


ra=l+efn 




E 

n=3fc+5+2^+£^. 


+ E E 

A :=0 ^n=3fc+3+^+£‘^^ 

Cd log'=+i(2/i?) 


eulerlogd^(^) 


R” 


- E 

71 —3A;-1-4-1-'^+£^^ 


B^^ eulerlog^(i?) 

^n+l /2 


i?" 


(25) 


[Recall that is defined in Eq. (1^01) .] In particular, note 
that the individual inodes are purely rational integer- 
order PN series until the first appearance of the loga¬ 
rithm, where they start to have transcendental contri¬ 
butions, as well. The transcendentals and the logarithm 
first only appear together in the form of eulerlog^(i?) 
in the integer-order terms, but then, starting with the 
\!R^+^+^i-'^ term (i.e., the same order at which eulerlog^ 
terms start to appear), they also get tt^ and C(3) terms, 
where C is the Riemann zeta function. Much of the in¬ 
crease of complexity is described by the simplifications 
[Eqs. (fT51) and (l^ j. though there are a few logarithms, 
transcendentals,^ and half-integer terms that the sim¬ 
plifications do not remove, just as found for the en¬ 
ergy flux in [ 3 ^, which should be expected, from the 
form of the expressions used to obtain both quantities. 
In particular, the expression in Eq. (l25|l does not in¬ 
clude the appearance of the eulerlog 2 (i?) log(2/i?)/i?^^ 
and TT log(2/i?)/i?^'^-® terms that are known in the (2,2) 
mode of At/ (which are also not given by either of the 
simplifications). 

We apply the PSLQ integer relation algorithm 
in its implementation as the FindIntegerNullVector func¬ 
tion in Mathematica to obtain analytic forms for the 
PN coefficients of the retarded (t, m) modes of At/.® 
Specifically, if one inputs a vector of decimals to the 
PSLQ algorithm, it returns the (nonzero) vector that is 
orthogonal to the input vector and has a small (L^) norm. 
One can thus apply PSLQ to identify the analytic form 
of a number from a sufficiently accurate decimal repre- 


Note that 7 and Q evaluated at odd integers are not known to be 
transcendental, or in most cases even irrational. However, they 
are all strongly conjectured to be transcendental, so we shall refer 
to them as such. 

® Note that we shall often use the name PSLQ as a shorthand for 
the FindIntegerNullVector function. 


sentation, if one knows (or has an educated guess for) the 
transcendental numbers [here, for instance, tt, log(2), 7, 
^(3), etc.[ present in the analytic form. This is particu¬ 
larly simple when the number one will obtain is a linear 
combination of the transcendentals with rational coeffi¬ 
cients, as is the case here, where the vector in question is 
simply the decimal expansion of the number to be iden¬ 
tified, along with 1, and any transcendentals thought to 
be present. Of course, PSLQ will give an output for any 
vector, but the outputs that do not correspond to a true 
relation are almost always large and “ugly-looking” for a 
sufficient number of digits, while the true vector will have 
a certain “nice-looking” structure (which we will discuss 
further later). 

C. An example: Obtaining the analytic form of the 
P-j coefficient of the full At/ from the decimal form 
given in SFW 

As an example, we consider obtaining the /Sy coeffi¬ 
cient [i.e., the coefficient of the log(i?)/i?® term in AU] 
from the numerical value given in SFW, as was reported 
there (and confirmed by the analytic calculation in [^1. 
This is the coefficient of log(i?) at the first order where 
there is a log^(i?) term, so, by analogy with the transcen¬ 
dentals appearing in the nonlogarithmic term at the first 
appearance of log(i?), we expect to have 7 and log(2) 
terms here. As we saw in the expression of the structure 
of the PN expansion of the modes in terms of eulerlogs 
above, this linking of log(i?), log(2), and 7 is a generic 
feature, though it is broken at high orders by the ap¬ 
pearance of the log(2/i?) terms. Indeed, since this is the 
first appearance of log^(i?), and thus only comes from 
the (2, 2) mode, we can actually subtract off the log(2) 
and 7 terms and only have to obtain the rational piece 
using PSLQ. However, we shall first consider the case of 
using PSLQ to obtain the full term, since this is how we 
initially obtained it. 
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Starting from 

Pr = 536.4052124710242868717895394750389112702062 
69552321207927883360240368736326766131833.. 

(26) 


which is taken directly from Table 1 in SFW, we can 
apply PSLQ in the form of Mathematica’s Findlnte- 
gerNullVector function to the vector {/Sy, 1, 7 , log(2)} and 
obtain the expression of 


„ 5163722519 109568 219136, 

Br = - 7 -log( 2 ) 

^ 5457375 525 ^ 525 ’ 


(27) 


with at least 42 digits. {Nota bene: We find that the 
final digit given by SFW is incorrect and should be a 
6 .) We were able to reject the expressions produced by 
smaller numbers of digits since they lead to anomalously 
large prime factors in the denominator (i.e., the term in 
the vector returned by PSLQ corresponding to ^7 itself), 
except if one only evaluates to a very small number 
of digits, of course: See Fig. [T] (cf. Fig. 5 in [s^, which 
shows an alternative method for detecting a likely true 
relation using PSLQ by looking at the size of the smallest 
entry in the vector versus the number of iterations of 
the algorithm, which is not information available when 
using Mathematica’s FindlntegerNullVector). We shall 
discuss this method of looking at the prime factorization 
further in Sec. EIB 



Number of digits in 


FIG. 1: The largest prime in the denominator of the expres¬ 
sion for Pj returned by Mathematica’s FindlntegerNullVec¬ 
tor function applied to the vector {Pt, 1, 7 , log(2)} for Pr eval¬ 
uated to varying numbers of digits from 10 to 45. The vertical 
dotted line marks the point at which this function returns an 
accurate analytic expression for P-j. 


Interestingly enough, the minimum number of digits 
required to obtain this expression accurately with Find¬ 
lntegerNullVector is somewhat dependent on the order of 
the terms in the vector. For instance, if we consider in¬ 
stead the vector {/ 37 ,log( 2 ), 1 , 7 }, we only need 39 digits 
to obtain an accurate analytic form. If we scale by the 


denominator of /Sg (which is 575), i.e., consider the vector 
{ 575 / 37 , l; 7 i log(2)} we also only need 39 digits; this sort 
of scaling is much more effective at higher orders where 
the denominators are much larger, and can decrease the 
minimum number of digits required to obtain an accu¬ 
rate expression by 27 digits or more. If we subtract off 
the 7 and log( 2 ) terms using the expectations from the 
eulerlog 2 (i?) form of the log^(i?) coefficient at this order, 
then we only need 22 digits to obtain the resulting Pt 
(here we use the vector {, 37 ,1}, of course): While Find¬ 
lntegerNullVector returns the correct result with between 
12 and 15 digits in this case, it then returns an erroneous 
relation when one uses between 16 and 21 digits before re¬ 
turning to the correct result for 22 digits and above. We 
have not observed such unusual behavior in our other 
determinations. Here scaling by 575 actually increases 
the minimum number of digits required for an accurate 
result by 1 , which can sometimes be the case when one 
is scaling by a relatively small number, as here. 


D. Combining together the values at different radii 
to increase the number of digits known 


For the low-order coefficients, which are purely ratio¬ 
nal integer-order PN coefficients, we can apply PSLQ di¬ 
rectly to an appropriate number of digits at a given radius 
(e.g., for R = 10®°, one can expect to get at least ^ 40 
accurate digits at a given order for the leading term). 
If one can identify the rational represented using PSLQ 
with this number of digits, then one can subtract it off 
and proceed to the next order. For the higher-.^ modes, 
the purely rational coefficients persist to high enough or¬ 
ders and are large enough that one needs more than the 
number of digits provided by merely evaluating A[/ at 
R = 10™, the largest radius we consider. In such cases, 
and also when we need to consider cases with logarithms 
and transcendentals at higher orders, we combine to¬ 
gether the values from several radii (up to as many as 
15 radii for certain high-order pieces). The expressions 
we use for this purpose are long and unilluminating, but 
we give a simple example here to illustrate the method. 

Let us assume that we are at a point in the computa¬ 
tion where we expect that the first few terms of the PN 
expansion of the mode we are considering to look like 


o lr>\ “-^.0 , Q^IV-l-1,0 + aA-l-1,1 log(i?) 

Sn{R) = ^ + -- 

+ 0{R-^-^) 


(y-N+2,0 

rN+2 

(28) 


(taking the R~^~‘^ term to have no logarithms, since 
we are just interested in its overall scaling, even though 
this will never be the case in this sort of situation in 
actuality), and we wish to obtain to ~ 2k digits. 
Now, since we know the value of Sn{R) at R = 10^’ for a 
range of integers fc, we can combine together the values 
of SNiR) at the radii R = 10'=, 10'=+^, and 10'=+« (with 
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p,q G N such that 10 ^+^* and 10 ^+'^ give radii at which we know the value of Sn), giving 

I 


aN,o 


{q - p)S'Ar(10'=) - gl0(^+i)P5Ar(10'=+P) + pl0(^+i)«S'jv(10'=+«) 
q — p — qlQP +plO'J 


q-p-qlO P+plO aN+ 2,0 
q-p-qlQP+pm 102 '= ■ ^ ’ 

' -V-' 

n 


Here the remainder TZ gives a small enough correction 
that the first term will give aNp to ~ 2k digits, pro¬ 
vided that aAr+ 2,0 is not much larger than aNp, which 
will usually be the case. (Note that we have neglected 
the rest of the remainder, whose leading term goes as 
10“^'=. We also have not included half-integer terms in 
the remainders for simplicity, though the presence of a 
half-integer term before the given remainder term will, 
of course, reduce the number of digits one obtains from 
the expression.) 

We obtained the particular linear combination given 
in Eq. (1^ by considering 5'jv(10'=) -I- ASn{W^~^p) + 
i35Ar(10'=+'^) and fixing the coefficients A and B by de¬ 
manding that the resulting expression does not contain 
the two R~^~^ terms (viz., ajv+i^o and aAr+iq). One 
then solves the resulting expression for to obtain 

Eq. (|29p . One derives the more involved expressions for 
more complicated cases with more terms and more radii 
in the same way by solving a linear system, which Math- 
EMATICA will do quite efficiently. For the determinations 
reported in this work, we needed at most 1240 digits and 
15 radii, which we used to determine the coefficient of 
the {R )/term in the (2,2) mode [and implicitly 
check the prediction of the simplification for the coeffi¬ 
cient of the log^(i?)/i?^^'^ term]; removing these terms 
was necessary to obtain the nonlogarithmic part of 12PN 
coefficient of the ( 2 , 2 ) mode. 


E. Overview of our method for obtaining analytic 
forms of the coefficients of the modes of AU 

Our general approach for determining analytic forms of 
the PN coefficients of the modes of AU is first to obtain 
the coefficient of the highest power of log(i?) present at a 
given order, which will always be rational (or a rational 
times TT, for the half-integer terms), and will always come 
from the corresponding power of eulerlog^(i?), where m 
is the mode’s degree (i.e., its magnetic quantum num¬ 
ber). Indeed, once we have obtained the coefficients of 
the first three log(i?) terms in the PN expansion of a 
given mode, we are able to predict the coefficients of all 
of the highest powers of log(i?) (and, in fact, much more) 
using the simplification. Thus, while an individual PN 
coefficient of the ( 2 , 2 ) mode of AU can have as many as 
17 transcendentals at the relatively high PN orders we are 
considering, we have to use at most 5 transcendentals in 


the vector to which we apply PSLQ (for the lOPN non¬ 
logarithmic term), since the coefficients of the remaining 
transcendentals are predicted by the simplification, or 
given by the coefficients of higher powers of logarithms 
at that order. (We only need 4 transcendentals for the 
nonlogarithmic piece at 12PN, despite its more compli¬ 
cated structure, since at this point in the calculation we 
have removed some of the transcendentals that entered 
at lOPN using the simplification.) 

Once we have obtained (or—more often—checked the 
simplification’s prediction for) the coefficient of the high¬ 
est power of log(i?) at a given PN order, we then 
subtract off the appropriate rational times a power of 
eulerlog^ (i?) and proceed to the lower powers of log(i?), 
which have a more complicated structure. At the orders 
where there are powers of log( 2 /i?) present, we still sub¬ 
tract off the putative contribution as if the log" (i?) term 
came solely from a eulerlog]]^ (i?) term, and then include 
the appropriate piece when obtaining the coefficient of 
the next lower power of log(i?) to account for the pres¬ 
ence of the log"(i?) term. For instance, when the log(i?) 
term comes from aeulerlog^(i?) -I- 51og(2/i?), so we sub¬ 
tract off {a + 2b) eulerlog^(i?), taking the log(i?) term to 
come solely from an eulerlog^ (i?) term, we thus include 
the remaining transcendental, viz., 27 -f log( 2 )-|- 21 og(m), 
in the vector to which we apply PSLQ when obtaining the 
coefficient of the nonlogarithmic term at this PN order. 

The only exception to this procedure occurs when we 
can predict the coefficient of the eulerlog]];^)!?) contri¬ 
bution from the simplification (necessarily for n > 2 , 
since the coefficients of the eulerlog^(i?) terms are in¬ 
puts to the simplification, and thus not predicted by 
it), in which case we simply obtain the coefficient of 
log"( 2 /i?) directly from the coefficient of log"(i?). Ad¬ 
ditionally, at 12PN in the (2, 2) mode, things are quite 
complicated, since we have to disentangle contributions 
from eulerlog 2 (i?), eulerlog 2 (i?), eulerlog 2 (i?) log( 2 /i?), 
log(2/i?), and log^(2/i?) terms. Here we just subtract 
off the log^(i?) coefficient we obtained as a log^(i?) term 
and then include 7 -I- log(2) and log(2) in the vector to 
which we apply PSLQ to obtain the log(i?) coefficient. 
The coefficients of 7 -I- log(2) and log(2) in the log(i?) 
coefficient then let us predict the coefficients of certain 
7 log( 2 ) and log^( 2 ) contributions in the nonlogarithmic 
coefficient, so we need to include only 27 - 1 - 31og(2) and 
72 -I- 37 log( 2 ) -|- (9/4)log^(2) in the vector to which we 
apply PSLQ. 
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For the more complicated terms, we use the modes of 
the energy flux at infinity as a guide to the transcenden- 
tals we expect to be present that are not already pre¬ 
dicted by the form of the simplification we have deter¬ 
mined at a given order. (These modes have been calcu¬ 
lated to 22PN by Fujita [^, with simplified and factor¬ 
ized forms given in [38|.) Once we have obtained a new 
transcendental (or other new contribution that looks as 
if it is part of the simplihcation) at a given order® in 
a few modes, we conjecture its general appearance in 
the simplification and then check it using a few other 
modes. Again, comparison with the analogous Sim fac¬ 
torization of the modes of the energy flux from [s^ is 
useful in determining what likely comes from the simpli¬ 
fication, though the structures are not exactly the same. 
We also may need to obtain simpler higher-order terms 
(i.e., higher powers of logarithms at higher PN orders, or 
half-integer terms, which are simpler than integer-order 
terms at comparable orders) first, in order to remove suf¬ 
ficient terms from the series so that we can obtain the ex¬ 
pression to enough digits with a relatively small number 
of radii. 

Finally, as mentioned above, we generally scale by the 
denominator of the rational term at the previous PN or¬ 
der (i.e., the current PN order minus 1), since we find 
that this significantly reduces the number of digits re¬ 
quired to make an accurate determination of the analytic 
form. One exception to this is any case where the simpli¬ 
fication predicts most of the term [e.g., up to a log"(2/i?) 
contribution, or the additional contributions to the non- 
logarithmic part of half-integer terms at higher orders]. 
In this case, we will not scale at all, or scale by the de¬ 
nominator of the similar addition at the previous PN 
order, as these additional terms have significantly sim¬ 
pler denominators than the full coefficient. With these 
scalings, we need at most 249 digits to successfully ob¬ 
tain an analytic form. This maximum is needed for the 
nonlogarithmic term of the (2,2) mode, the most 

complicated coefficient of a mode we consider. 

For the modes with £> 4, we also scale by a high 
power of m (as high as 12 for ^ = 9 and 10) when deter¬ 
mining the (linear) log(i?) and half-integer terms, since 
such high powers of m occur there. For these modes, we 
are able to use the simplification to make the determi¬ 
nation of the analytic form of the PN coefficients mostly 
automatic (in the sense that we only have to choose the 
number of digits used and various scalings, and verify 
that the results satisfy all the consistency checks we dis¬ 
cuss below). Note also that we are only considering the 
linear log(i?) terms here, since all the higher powers of 
log(i?) in these modes are given by the simplification, to 


Here we consider not the absolute PN order, but the relative 
PN order past the first appearance of an eulerlog term, since a 
term at such a relative PN order will have the same complexity 
in all modes, while the complexity at a given absolute PN order 
decreases with increasing i + eim ■ 


the order we are currently working. 


F. Checks on the output of PSLQ 

When we are performing these PSLQ determinations, 
it is important to ensure that one has sufficient accuracy 
(and the correct transcendentals in the vector) so that the 
form returned by PSLQ is reliable. Besides the basic test 
of making sure that the analytic form does not change 
as one increases the number of digits, up to the maxi¬ 
mum number that are expected to be given accurately 
by the combination of radii being used, a very stringent 
test is generally making sure that the denominators do 
not have any large prime factors. Such smooth numbers 
are distributed relatively sparsely among large numbers 
when the largest prime allowed is relatively small (e.g., 
smaller than the logarithm of the large number): See, 
e.g., Granville’s review for a discussion of the prop¬ 
erties of smooth numbers. 

In particular, Granville gives an upper bound on the 
smooth number counting function in his Eq. (1.23) that 
gives an easy way to see how unlikely it is for a randomly 
selected d digit number to have all its prime factors less 
than p. This probability will be less than 

1 fljd + 1) log(lO)/ log(2)J -b 7r(p) 

IQd-fi _ IQd 

where (’) is the binomial coefficient, [•] is the floor func¬ 
tion, and 7r(p) is the number of primes < p. This prob¬ 
ability is generally extremely small for the cases in ques¬ 
tion. For instance, the denominators of the purely ratio¬ 
nal terms at 12 and 12.5PN each have 32 digits, but their 
largest primes are both 19, for which the probability is 
less than 10“^^, for a randomly selected 32 digit number. 

In addition to making sure that the denominator is a 
smooth number, other consistency tests include check¬ 
ing that the result is insensitive to small changes in the 
prime factorization of the overall scaling, or computing 
the coefficients of (7“A[/ for different a (e.g., a = —1/2 
and a = -1-1/2, though we used other values, as well) and 
making sure that the results are consistent. One also ex¬ 
pects that terms that have only recently started to appear 
in the expansion will have simpler forms than those that 
have been present in the expansion for longer, which also 
allows one to reject some spurious expressions. In partic¬ 
ular, one expects to see powers of the characteristic large 
prime from the numerator of the first PN coefficient of v 
(cf. Table 1 in [^) in such terms that have just started 
to appear. 

Finally, if one happens to have many digits for a given 
term, other consistency checks include adding on other 
transcendentals to the vector and making sure that PSLQ 
gives zero coefficients for them, or obtaining the result 
without subtracting off some of the known results from 
the simplification (particularly if this is a term one al¬ 
ready has to include in the vector anyway, since the sim¬ 
plification only gives a portion of it). All these checks, 
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plus the overall consistency check that we continue to 
obtain expressions of the expected form for the vari¬ 
ous modes (given the simplification and the form of the 
modes of the energy flux) as we continue to high orders, 
give us high confidence that our analytic results for the 
modes (and, as described later, the full AC/) are indeed 
the true ones. This confidence was recently confirmed by 
the exact agreement of our results with the 20.5PN re¬ 
sults of Kavanagh, Ottewill, and Warded [l^, obtained 
completely analytically. 

IV. TERMS PREDICTED BY THE 
SIMPLIFYING FACTORIZATION TO 
ARBITRARILY HIGH ORDERS 

The two simplifications we have found also predict cer¬ 
tain higher-order logarithmic and half-integer terms, ex¬ 
tending to arbitrarily high orders, since we assume that 


the and portions of the 

simplifications hold to all orders, as we expect them to, 
since the similar Stm and 14 m factorizations found for the 
energy flux in [s^ hold to high orders (presumably to all 
orders) and contain the same exponential terms. More¬ 
over, we can see where these factors arise in the calcula¬ 
tion from a study of the structure of the MST solution 
used to compute AC/; see the Appendix. In particular, we 
can predict the coefficients of the first five appearances 
of a given power of log(i?) in both the integer-order and 
half-integer terms. The only thing that prevents us from 
being able to predict further terms is the appearance of 
pieces that are not given by these simplifications at higher 
orders. 


Specifically, the higher-order logarithmic and half¬ 
integer terms in the full AC/ that the simplification pre¬ 
dicts are given by the appropriate terms from 
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+ {18PN and higher leading logarithmic terms} + {other terms that are not accurate predictions}, 
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3 
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(32b) 


and the “appropriate terms” that one should take from 
this expression are (as discussed above) the first five ap¬ 
pearances of a given power of the logarithm in each of 
the integer and half-integer PN coefficients. (Note that 
the final subtracted terms in Cf^ and Cfj^ are only nec¬ 
essary to remove some low-order nonlogarithinic integer- 
order terms, so one could leave off the subtracted term 
and simply ignore the additional terms, since they do 
not mix with the predictions of the simplihcation.) While 
this expression gives further appearances of these powers, 
such further appearances are not predictions for complete 


J 


coefficients of the full At/, which obtains contributions 
from terms that are not included in the simplification. 
(We give code that picks out the accurate predictions in 
. the Supplemental Material [i^.) Moreover, even if there 
were no additional terms besides those given by the sim¬ 
plification, one would need to add on more terms in the 
series for a given mode, as well as additional modes, to 
obtain the sixth and higher appearances of a given power 
of a logarithm. Also, note that here we only need f/^m 
to second order [i.e., to 0(1/77®); recall that u is a series 
in 1 / 778 ] in Cf^ and to first order [i.e., to 0(1/773)] in 
Cfm I ^’^d these higher-order terms are only needed for the 
highest few appearances of a given power of a logarithm. 


It is also possible to obtain explicit expressions for the 
coefficients of the first few appearances of a given power 
of log(77) in the PN expansion of At/ from Eq. (IHTl) . us¬ 
ing the Taylor series for the exponential function. As an 
illustration, we give here the expressions for the coeffi¬ 
cients of the first two appearances of the nth power of 
log(77) in both the integer and half-integer PN terms of 
At/: 


f64 

y856^ 

n—1 

1 , 





2816 /856y " 


16 /214 
45 vT05 


16 /856y+" TT 


704 /856y+" 8 /214 

TosVTosy ”45 1^105 


Here one must only consider n > 1 in the integer-order 
terms, but can consider n = 0 in the half-integer terms. 
Note that the first appearance of a given power of log(77) 
(in either the integer or half-integer PN terms) comes 
solely from the leading term of the series multiplying C 22 
(i.e., just from the dominant 2, 2 mode), while the second 
appearance comes from the next term in the series mul¬ 
tiplying C 22 in addition to the leading terms multiplying 
Cf^iov (7,m)e{(2,l),(3,3),(3,l)}. 
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243 /78\ 


243 /78y+" 
"84” 


1 /26Y 

1 /26y+" 
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773 


777-5 


log" (if) 
n\R?‘^ 

(33) 


One can also similarly predict arbitrarily high-order 
“leading logarithmic” terms in the energy flux using the 
St,rn factorization from [s^, though this was not noted 
there: Here one obtains the coefficients of the first six 
occurrences of each power of a logarithm in the integer- 
order PN terms and the first five occurrences in the half¬ 
integer PN terms. We give code that computes these pre¬ 
dictions in the Supplemental Material . Additionally, 
note that Nickel makes similar predictions of leading 
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logarithmic-type terms to arbitrarily high orders for the 
ground state energy of H 2 , and is able to derive some 
of them. Finally, these sorts of predictions of leading 
logarithms to arbitrarily high orders are likely related to 
the multipole moment beta functions discussed by Gold- 
berger et al. [H, , where they predict the coefficient 

of the first occurrence of a given power of a logarithm in 
both the energy flux and the binding energy using the 
beta function for the dominant (2, 2) mode. 


V. COMPUTING THE INFINITE SUM OVER 
RENORMALIZED £ MODES TO OBTAIN THE 
FINAL RESULT FOR AU 


AC/, being a conservative invariant (and thus coming 
from the half-retarded plus half-advanced field, as dis¬ 
cussed in Sec. 5 of @), requires a renormalization pro¬ 
cedure where a noncontributing singular part of the re¬ 
tarded AU calculated at the position of the particle needs 
to be subtracted. This is done by using a mode-sum reg¬ 
ularization technique where the retarded part is written 
as a sum over angular harmonics and the singular part is 
written as a sum over angular harmonics by extending it 
on the coordinate 2-sphere passing through the particle 
(i.e., at r = tq). The explicit equations used in renormal¬ 
ization are given in . The sum over i modes converges 
quite slowly (the summand goes as so it is custom¬ 
ary to improve the convergence by finding higher-order 
regularization coefficients numerically, as in [3, [52|, [6^ . 
However, to obtain an accuracy of N digits in the fi¬ 
nal result of the renormalized AU, one has to obtain 
these higher-order regularization coefficients to N digits 
as well, which would necessitate going to prohibitively 
large £ (e.g., /max 10^ for an accuracy of 5000 digits). 
(The Wentzel-Kramers-Brillouin method detailed in 
could be useful here in future work.) 

Nevertheless, it is possible to obtain the analytic form 
for a given nonlogarithmic integer order PN term by ob¬ 
taining the general form of the PN coefficient of a renor¬ 
malized I mode and performing the sum analytically, as 
in Bini and Damour [23 - [^ . We thus note that the 
nPN coefficients of all the renormalized I modes with 
/ > n — 1 are purely rational (with no transcendentals 
or logarithms). (Here we only consider n £ N, since the 
half-integer terms only have a finite number of I modes 
contributing.) One can thus easily obtain analytic forms 
for these coefficients using PSLQ. (Here one scales with 
the denominator of the previous I mode, to help with the 
determination, and needs at most the values at four radii 
to obtain enough digits at 12PN.) We then find that the 
general form of these coefficients as a function of I can be 
expressed as linear combinations of members of a small 


family of functions, namely 

(34a) 
(34b) 
(34c) 

[Note that V" is the only one of these where the effect 
of the superscript n is the same as taking V to the nth 
power.] These functions are similar to, though slightly 
more complicated than, the form considered for the reg¬ 
ularization coefficients in Sec. V of Shah et al. [1^. We 
solve the linear system to obtain the coefficients (noting 
that one obtains excessively large rationale for the coeffi¬ 
cients if one does not include the correct functions in the 
solve) and check that the expression successfully repro¬ 
duces the values of the coefficients that were not used in 
the solve. We need to go to / = 87 (starting from / = 16, 
to avoid logarithmic terms at higher orders) for 12PN, 
the most complicated case we consider, for a total of 72 
/-modes. 

The general expressions for the first six PN coefficients 
have the form 


r,"(/) := 
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V"(/) := 
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-b 


(-ir 


{i + k+ 1/2)^ {£ - k + 1/2)^’ 

1 , (- 1 )" 

(/-bfc)" (i-k + l)^’ 
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(/-b 1/2)"' 
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'1 5 

ri^_4 & Ul_s & V^, (35) 

& r^_2 & & ul_^ & Ul & Ul & V^, 

& T[^_3 & kTi k U\_^ k Ul_2 k Ul & 

where we just give the functions present, not the coeffi¬ 
cients, and a range in a subscript indicates that all the 
functions in that range are present. We give the explicit 
expressions up to 12PN in the electronic Supplemental 
Material [6^ , and only note here that the specifics of the 
functions present grows in about the way one would ex¬ 
pect: At nPN, one has Wf_(„_ 3 fe+ 2 ) te™s 

present, for fcs with n — 3k + 3 > 1 and n — 3fc -b 2 > 1, 
respectively, as well as Tp , Up terms for larger p (where 
the specifics of the terms present at a given PN order 
has a somewhat more complicated structure). One also 
has present for all even k < {2/3)n. For instance, at 
12PN, we have 

Ti^-i 2 & Ti "-9 & 7 )-6 & 7 )-3 kr^ k Ti® k rj k 
k U\_^^ k Ul_^ k U\_^ k Ul_2 k U\ k Uf k Ul 
& & V® & V®. 

(36) 

Also note that these general expressions diverge at the 
/s for which the PN coefficient is no longer purely ra¬ 
tional, due to the Ul^_i term (cf. the discussion of the 
appearance of the logarithms at places where there are 
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apparent poles in I in the general form in Sec. II F of 
Bini and Damour (^1. 

As Bini and Damour mention [^ . the infinite sums 
over these functions are straightforward to evaluate if 
one makes a partial fraction decomposition (and Mathe- 
MATICA will do them automatically without even needing 
to perform a partial fraction decomposition first): One 
finds that the sums over many of the terms telescope to 
a finite sum and the rest can be evaluated in terms of the 
Riemann zeta function evaluated at even integers (giving 
even powers of tt with rational coefficients). Since these 
general expressions are not valid for the low-^ modes, 
where there are also transcendentals present, one adds on 
the contributions from these low-order modes separately 
to obtain the final expression. One finds that the size of 
the numerator and denominator of the final purely ratio¬ 
nal term is a good indicator of errors in the calculation: 
If one has omitted a piece, or determined its analytic 
form incorrectly, this rational will be more complex than 
one would expect it to be, given the complexity at the 
previous order. 


While we have analytic forms of the PN coefficients 
for AC/ through 12.5PN, with the 13.5PN term and all 
but the nonlogarithinic piece of the 13PN also known, 
we only give the full AC/ to 11.5PN here to save space. 
The analytic forms of these high-order coefficients are 
quite lengthy, even when written in eulerlog form. We 
give the previously known lower orders (which we have 
re-obtained) in their eulerlog form as well, for compar¬ 
ison, and to illustrate the structure. We also give the 
expression with the terms given by the simplifications 
removed, where we go all the way to 12.5PN. We give 
the full expressions for all these quantities in the elec¬ 
tronic Supplemental Material [b^. Here we scale AC/ by 
u := 1 /i? and write the expansion in terms of u, so that 
the coefficient of u" gives the nPN term of AC/. We also 
abuse notation (i.e., we use “physicist’s function defini¬ 
tions,” not “mathematician’s function definitions”) and 
write eulerlog,„('u) := 7 -t- log( 2 mu^/^), which has the 
same value as the previous expression in terms of R if 
one uses the u related to this R, but is, of course, not 
given by substituting u for R in the previous expression. 
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24640 
82561159 
467775 


-TTU 


7.5 


31185 
7516581717416867 
34763478750 


19008 
246847155756529 


4 41408,64, , 31988738821 , , , , 

' - - y - 1222452000 <"> 


18496880640 

80813099648 


eulerlog 2 (u) 


6442450944 105 ' 5 ^ 1222452000 oiw ' 33108075 

85126268709 , , , , 67792273408 , , , , 798828125 , , , , 3359232 , , , , 

eulerlog 3 (M)-eulerlog 4 (M) -k eulerlog 5 (u) - eulerlog 6 (u) 


15695680 
16022 


eulerlog! (it) — 


70945875 

4820992 , . 2 . . 18954 


11025 
23447552 
55125 
32962327798317273 


eulerlog 2 (it) 


11025 


, 8.5 


■ eulerlog 2 (u) -k 


49 


741312 
eulerlog 3 (it) 


17875 
2207224641326123 


10480362137370508214933 


1048863816000 
11665762236240841 


219136 
1575 ’ 


2044301131372500 


226072985600 


TT^ — 


27101981341 


-7r° -k 


10221088 


448 


C(3)-— log(2u) - 


61470271483 


eulerlog! (u) 


549755813888 100663296 2835 ^ 5 ' 814968000 

2840603616267776 , , , , 8677864251603 , , , , 5946112890241024 , , , , 

eulerlog2(u) H- onriononnn — eulerlog3(u)- AAr,Aor\Anno'7r : — eulerlog4(u) 


442489422375 
58533203125 


15567552 

5373212672 


eulerlog5(M) 


392392000 
309049344 


9823275 
1055996 , 


eulerlog 2 (u) — 


49 


125125 

eulerlog3('u) 


eulerlogg(u) — 


442489422375 

96889010407 , , , , 51178 , , 2 / ^ 

eulerlog 7 (u) -k , eulerlog!(u) 


277992000 


2 . , 94770 , , 2 , ! 1647312896 , , 2 


1964655 


eulerlog! (u) it® -k 


19845 
30185191523470507 


12236744520000 


11025 


1712534 

1157625 


eulerlog! (it) — 


1031692288 

1157625 


eulerlog 2 (it) 


246402 

343 


eulerlog 3 (it) 


, 9.5 


















































































18 


238946786344653264799175280203 1070208441923650860489683 


4522423405833558225000 


58656715985387520000 


832229033014028790267991_ 
1662461581197312000 ’ 


+ 


+ 


54067065388369 g 128695611256 

•TT — 


12884901888 ' 5457375 ^ 

1483437716511288604288 46895104 


32768,,,, 20416, , 157982464536376957 , , , , 

-c -oc— -eulerlogi(M) 


35 


674943865596000 


14480466347221875 
8040008069311889408 


33075 


3506176 
525 


C(3) eulerlog 2 (u) - 


52813127885844357 

10492954472000 


eulerlogg (u) 


, , , , 263296063591796875 , , , , 13640920722432 , , , , 

eulerlog 4 (w)- o^.oi eulerlogg (m) -eulerlogg (u) 


82745521984125 ^ 8742130068672 

6491563697269 , , , , 1099511627776 , , , , 

eulerlog 7 (u)- i.no.noor; eulerlogg (u) 


1181466000 
69907855522816 
3781960875 
7548828125 


eulerlogg (u) 


1688511825 

79338802833 

61661600 


1146520375 
eulerlogg (u) 


eulerlogg (m) — 


+ 


4077216 

1096738 


eulerlogg (m) — 


187580416 

165375 


eulerlogg (u) 


,10 


+ 


54944178599 
' 7491884400 
705049919488 , , 2 . ^ 

108056025 

54441085537326639211 


3824681905044000 


+ 


78847804_ 
66825 ^ 


416745 
10351714238464 
6807529575 


, , , , 38224121176064 , , , , 1232010 , , , , 

eulerlogg (m) - o/mo-yd-yo-yc eulerlog2 (u)-— eulerlogg (u) 


34037647875 


eulerlog 4 (it) 


,10.5 


343 

3814229145040080910470246242071097 


13798711885916862654750000 
1213451006696869077146724173 


497508986166915487823810257447 2 

36601790774881812480000 ^ ^ 2234348365129187328000 

27594682424243283328.,., 175005952, , 54784, 2 


— 


5276940898567193189_ 

25975962206208 


-C(3) 


C(5)- 


+ 

+ 


19864845 ^ 105 

292720019838735815778069367 


55125 
1712534 


313683671243842099200000 

694575 ' 

17152255889408499680050304 

2063384576 

325477442086506084375 

694575 

37422973611649093363871733 

1232010 . 

-1--TT' 


■ log( 2 u) 
128176 


525 


^2- 


11025 
154271744 


11025 
1364688 


log"( 2 u) 

C(3)j eulerlogg (u) 

C(3) ) eulerlog 2 (u) 


1595162202161082218708992 

9299355488185888125 
12192267599501090688 


eulerlog 4 (u) 


49 


282725878023723294921875 


202481230826875 
690493302243328 


eulerlogg (u) — 


829173552753401856 
60335405593501190803 


C(3) ) eulerlogg (it) 

eulerlogg (u) 


57747104415 
501413283015334912 


eulerlogg (u) — 


1792954994688000 
205891132094649 


eulerlogg (m) 


168551219200 
12234781198165473 

196392196000 


eulerlogg (u) 


eulerlogg (u) 

792734736884113 


14316991088400 


eulerlogg (u) 


eulerlogg (u) 


1768982113681408 

127830277575 


eulerlogg (u) 


22370298575625 

3087470703125 , , ^ 72640032768 , , x 6850136 , , g, , 8253538304 , , g, , 

■ eulerlogg(u) + eulerlogg(u) - eulerlogg(u) + eulerlog2(u) 


+ 

+ 


159011424 
985608 , , 3 , ; 

3^3 eulerlogg (u) 

375160832 


3472875 


3472875 


45399846479271440442297518687 10107325522351333 


663973981856472412125000 


1311079770000 


-TT^ + 


3506176 
23625 ’ 


55125 
732046976712531 


839591622096533 , , , , /2629370415206008832 375160832 , 1 , , / n 

C(3) H-eulerlogg (it) +- 1 -tt eulerlogn(it) 

’ 112490644266000 ’ \ 65522472159375 165375 ' ’ 


308616308000 
20071104512 


, , , , 4430533694062592 , , , , 5835244140625 , , , , 

eulerlogg (it)-eulerlogg (it) H-eulerlogg (u) 


374414126625 


1749125664 


5788125 


eulerlog 2 (it) 


7rit^^° + 0{u^^). 


(37) 


The terms through 13.5PN for which we obtained analytic forms that we do not show here (i.e., all of these terms 
except for the nonlogarithmic 13PN term) have the expected increase in complexity, given the pattern at lower orders, 
and the complexity of the energy flux at infinity (see |38l. l39l|). In particular, we obtain a tt® term (from the sum over 
all ^ modes) at 12PN, along with an eulerlogg (it) log(2it) term [from the (2,2) mode alone]. We also see the expected 
TT^ and C(5) terms in the linear logarithmic term at 13PN (which we obtain from our fit to the full AC/, as described 
below) and get the first log(2u) term in a half-integer piece at 13.5PN. 
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If we write AC/ as a remainder plus the terms given by the two simplihcations, we have, now going all the way to 
12.5PN, 


AC/ 

u 


-1-2m-5u^ + 


121 


41 


1157 677 


512 


15 


1606877 60343 


3150 


768 


+ 


17083661 1246056911 


2800873 


4050 1769472 262144 

7516581717416867 741674600438227 


+ 


12624956532163 8826302018257 


382016250 


2477260800 


^2- 


23851025 

16777216^ 


22759807747673 


34763478750 


55490641920 


6442450944 


10480362137370508214933 25850880135908623907 


2044301131372500 


494421619507200 


-TT^ + 


32962327798317273 
549755813888 


— 


27101981341_ 
100663296 ’ 


+ 


238946786344653264799175280203 18927900985563784496607041 


4522423405833558225000 


1583731331605463040000 


166464310992954697947563 _ 
332492316239462400 ’ 


54067065388369_ 
12884901888 ^ 




8192 


-TTM 


10.5 


3814229145040080910470246242071097 31306151918920018820376206813081 


13798711885916862654750000 


2305912818817554186240000 


48548370032386602537826133 
89373934605167493120 


— 


5276940898567193189_ 
25975962206208 ^ 




2192896 
6741 


-TTU 


11.5 


+ 


176538264526096039025674251863176559931511 230595271714866856972896270561913344391 


12080725340499801136900598850000 


558842583466071892143636480000 


902366950567138330511081959572149 
18875774988611374546944000 


■k'^ — 


65536 

75 


log(2u) eulerlog 2 (u) 




31486592 

43335 


167517791710253563186933 

26599385299156992 
10 i 

12.5 I \ ' \ ' /^[l] "<"51 


-TT 

4 


-TTM 


E E + E E 


44336492264184971 
2473901162496 ’ 

[2] ^S1 


m—1 


i—2 m—1 


(38) 


where Tf^ and are integer order power series in 
u with rational coefficients, which we give (to the order 
known) in the electronic Supplemental Material [see 
Eqs. (|7T|l and (l24)l for the expressions for the (2,2) mode 
of AU/U]. (Note that the odd m terms for the C = 10 
Tf^s and the C = 4 Tf^s do not contribute until 13PN.) 

We find that the 13.5PN piece of AC/ has more terms 
that are not removed by the simplification than do the 
previous half-integer PN terms, just as occurs at this or¬ 
der in the energy flux (see the expression for the Sim 
factorisation of 1722 in the electronic Supplemental Mate¬ 
rial for [ 3 ^), and, as in the energy flux, the additional 
terms all come from the dominant (2, 2) mode at this or¬ 
der. Specifically, the 13.5PN piece of AC //u remaining 
after subtracting off the parts given by the simplification 
is 


2096793662144 131072 


-b 


139033125 
7012352 


2625 


log(2M) 


225 


,13.5 


14024704 

7875 


eulerlog 2 (u) 


(39) 


However, the portion remaining in other PN coefficients 
of AC/ after using the simplification does not have exactly 


the same structure as that in r722/|5'22p- For instance, 
i?22/|5'22p also has eulerlog 2 and eulerlog 2 terms in the 
12PN coefficient. 


A. Checking the results for AC/ by making an 
independent fit 

We performed an independent check of these results by 
making a fit for the PN coefficients of AC/ using data at 
smaller radii and the fit procedure described in SFW [2lj . 
In addition to checking the decimal expansions of the 
terms we have already obtained analytically, we also im¬ 
plicitly check all the coefficients we have obtained in the 
fit by verifying that the higher-order coefficients are not 
too large, as described below. We perform these fits it¬ 
eratively, obtaining analytic forms for as many terms as 
possible with the accuracy obtained from a given fit, sub¬ 
tracting these off, and fitting again. In this case, we 
proceeded through six iterations, where the first fit only 
went to 20PN, and these coefficients were obtained with 
an accuracy of just a few digits, while at the fifth and fi¬ 
nal iteration, after we had subtracted off 48 coefficients, 
we obtained the 20PN coefficients that we did not ob¬ 
tain analytically to ~ 41 digits, and were able to go all 
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the way to 21.5PN, where we obtained the coefficients we 
did not know analytically from the simplification to ~ 10 
digits. 

We made verifications of these results by checking that 
the analytic forms we obtain have the expected forms, 
and that the terms given by the simplification agree, in 
addition to the stringent verification provided by the fit 
itself, described below. We used the simplification to aid 
this process, so we needed to include at most 3 transcen- 
dentals in the vector to which we apply PSLQ (for the 
16.5PN linear logarithmic term). This procedure (of us¬ 
ing PSLQ to iteratively improve a fit, aided by a conjec¬ 
ture for the form of certain leading logarithm-type terms) 
is very similar to the one used by Nickel to obtain high- 
order terms in the expansion of the ground state energy 

ofiZ+inlig. 

We give the final results of this fit (both analytical and 
numerical) in the electronic Supplemental Material [^ . 
including showing the remainder of the analytic terms 
after removing the portions given by the simplification. 

One advantage of using high-precision data to extract 
PN coefficients is that it is relatively easy to check the 
accuracy of the analytical coefficients calculated using 
PSLQ. If we had used an incorrect coefficient, say for 
an nPN nonlogarithmic term, and used it to find other 
coefficients, the coefficient of the nPN higher logarithmic 
terms and subsequent higher order PN coefficients would 
have increased by many orders of magnitude to give a 
nonsensical result. 

Let us illustrate this with an example. The numeri¬ 
cally extracted 21.5PN nonlogarithmic coefficient (0:21.5, 
in the terminology of SFW) has a size of ~ 10^^. If we 
had used an incorrect 2IPN log®(i?) term (C21), the 021.5 
we extracted from the fit would have increased to a size 
of about 10"^°, a nonsensical result. The analytical form 
of <^21 was determined from its numerical value, which 
was extracted with an accuracy of 13 significant digits. 
So, to test the sensitivity of the fit to the values of the 
digits we did not extract, we inject random analytical 
(absolute) errors in (,21 of magnitude ranging from 10“^^ 
to 10“^”^ and extract 021.5. These errors are injected by 
using random numbers between 1000 and 5000, multi¬ 
plied with powers of 10 ranging from —16 to —30. We 
see that if we had included an error of magnitude 10“^^, 
the numerically extracted 021.5 would have had a size of 
^ 10^^, and if we had included an error of magnitude 
10“^^ (which is more than twice the number of signifi¬ 
cant digits used to calculate the analytical form of C21), 
021.5 would have had a size of ~ 10^®. 

This example clearly demonstrates the sensitivity of 
the numerical fitting technique we use and how the an¬ 
alytical forms of numerically extracted PN coefficients 
can be checked by injecting errors. Of course, it is al¬ 
ways possible to have a quantity that only differs from a 
reasonable-looking analytic form at extremely high posi¬ 
tions in its decimal expansion (see some of the examples 
given by Bailey and Borwein (d^ . Id^l. However, this 
seems quite unlikely to be the case here, particularly be¬ 


cause we have a good idea of the form of the coefficients 
and the growth of their complexity, from the forms of 
lower orders and the PN expansion of the energy flux at 
infinity. 

VI. CONVERGENCE 

It is interesting to consider the convergence of the high- 
order PN expression we have obtained for AC/. In Fig. 
we compare the convergence of the plain 21.5PN expan¬ 
sion of At/ with various resummations. Here we compare 
with the numerical data from Table HI in Dolan et al. 
for radii of 14, 5.6,10}M and with data from Table IX in 
Akcay et al. [ig for a radius of (I0/3)M ~ 3.33M, con¬ 
verting their h^,^(x) into our AC/ using their Eq. (17) 
and Eq. (2) in [^. We find that while the rate of conver¬ 
gence decreases as the radius of the orbit decreases (as 
expected), the series still converges reasonably well inside 
the innermost stable circular orbit (ISCO) at r = 6M, 
and continues to converge quite monotonically close to 
the fight ring at r = 3M, albeit extremely slowly. More¬ 
over, the exponential resummation (of the entire series, 
as originally proposed by Isoyama et al. [t^, not mode- 
by-mode, as in [3^ [7l|) improves the convergence sub¬ 
stantially for low to medium orders, particularly within 
the ISCO, though it makes it significantly less monotonic, 
and actually worsens the convergence at high orders in 
the strong field regime. 

If one performs a partial mode-by-mode exponential 
resummation, either exponentially resumming the modes 
through i = 10 and the remainder of the full AU sep¬ 
arately, or using the simplifications on the modes and 
exponentially resumming the portions that multiply the 
simplifications, as well as the remainders of the modes, 
this does not perform better than exponential resumma¬ 
tion applied to the entire expression (though it also does 
not behave as erratically as full exponential resummation 
at high orders in the strong field). If one just applies 
exponential resummation to the individual modes, then 
one finds that it does improve the convergence of some 
modes, particularly the ones with larger i — m. Factoring 
out the test particle binding energy, as done in Akcay et 
al. [ 1 ^ also improves the convergence, particularly near 
the light ring (where the test particle binding energy di¬ 
verges), but does not improve the convergence nearly as 
much as the exponential resummation on its own. 

One can also estimate the radius of convergence (in v) 
of the PN series for AU by looking at where a„ 

is the nonlogarithmic coefficient of u". If the series has 
no logarithmic terms, then the radius of convergence, Vr, 
is given by 1/vr = limsup„^oo al/'^. Thus, since AU di¬ 
verges at the light ring (as discussed in, e.g., 0) , one 
expects the radius of convergence of the PN series for AU 
one estimates in this way by taking only the nonlogarith¬ 
mic portion of the coefficients (and just considering the 
known orders, without taking the limit) to be close to the 
fight ring. Or, to put it another way, one expects the size 
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FIG. 2: Convergence of the 21.5PN expression for AU for orbits at various radii, comparing with the numerical data from 
Dolan et al. [3^ and Akcay et al. 0- Specifically, we show the convergence of the plain series, as well as the results of factoring 
out the test particle binding energy and/or performing exponential resummation on the entire series. 


of the coefhcients of u” to grow approximately like 3"^^. 
One indeed sees that this is the case, as shown in Fig. [3l 
where we plot the Schwarzschild radial coordinate of the 
radius of convergence estimated in this fashion. We also 
plot the growth of the PN coefhcients of the energy hux 
at inhnity, for comparison. 


VII. CONCLUSIONS AND OUTLOOK 

We have introduced a method for obtaining analytic 
forms of high-order post-Newtonian coefficients to lin¬ 
ear order in the mass ratio from high-accuracy numerical 
results from black hole perturbation theory. We have 
also given the hrst application of this method to the case 
of Detweiler’s redshift invariant, which (when evaluated 
in linear black hole perturbation theory) gives the lin- 
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FIG. 3: An estimate of the Schwarzschild radial coordinate 
(in M) of the radius of convergence of the PN series for AU/u, 
obtained from , where a„ denotes the nonlogarithmic co¬ 
efficient of v". We also show the same estimate for the test 
particle energy flux at infinity {dE/dt)oo, scaled by the New¬ 
tonian energy flux (from js^), for comparison. 


ear in mass ratio piece of the binary’s binding energy 
and the EOB radial potential. Here we have found ana¬ 
lytic forms for all these coefficients to 12.5PN, and have 
obtained mixed analytic-numerical results to 21.5PN (in¬ 
cluding analytic forms for the complete 13.5PN term, and 
all but the nonlogarithmic piece of the 13PN term), sub¬ 
stantially improving on the previous 9.5PN knowledge of 
this quantity. We also found a simplification of the indi¬ 
vidual modeSjSimilar to that found for the energy flux 
at infinity in 38|, which also allows us to predict certain 


leading logarithmic-type terms to all orders in the full 
AU. 


The new terms we have obtained improve the accuracy 
of the series, even inside the ISCO and near the light ring 
(though the convergence there is very slow, as expected); 
factoring out the energy, which diverges at the light ring, 
improves the convergence somewhat. Since exponential 
resummation of the individual modes of radiative quanti¬ 
ties improves the convergence much more than exponen¬ 
tial resummation of the full quantity (see 1^), we 
had hoped that there might be a better way of performing 
the exponential resummation here, which would behave 
better in the strong-held regime. However, our experi¬ 
ments in this regard were unsuccessful, in that we only 
obtained very modest improvements, much less than the 
best improvement of exponential resummation applied to 
the full series, though the improvements did not have the 
full exponential resummation’s erratic behavior. 

It might also be possible to use these high-order per¬ 
turbative results to improve convergence by Hnding non- 
perturbative pieces, using resurgence (see, e.g., [t^ for an 
application of these ideas in quantum mechanics). An¬ 
other possibility would be to try to resum the purely 
integer-order PN series with rational coefficients that en¬ 
ter into the simplihcation or its remainder, as was done 


for a (likely considerably simpler) self-force series in . 

We are now in a position to apply this method to the 
much more difficult case of perturbations of the Kerr met¬ 
ric. Here we will likely combine a study of AU with 
a study of the structure of the en ergy flux at infinity 
(computed numerically to 20PN in and analytically 
to IIPN in [7l|), since our previous study of this struc¬ 
ture in the Schwarzschild case was very useful in the 
present calculation. 
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Appendix: Obtaining the and 

g 2 vim iog(2/H) contributions to the simplifications of 
the modes of AU 

Just as one can see where the Sim and Vim factoriza¬ 
tions of the energy flux from arise from the MST 
formalism (as discussed in Sec. IV of [l^), it should be 
possible to see how the simplifications for AU we have 
found [Eqs. (flSll and (1251) ] arise from the MST formal¬ 
ism, and (in the best case) predict higher-order terms 
in them. However, we shall see that the situation for 
AU is more complicated than that for the energy flux, 
and will at present content ourselves with seeing how the 
g 2 p,,„euieriog„(K) iog( 2 /fl) contributions to the 

simplifications arise. Note that here we shall expand in v 
instead of R, for simplicity, and to avoid confusion with 
some other quantities named R. 

Specifically, if one looks at Eq. (29) in and our 
Eqs. m and dHI, one finds that the modes of AU have 
the form 




VE[i?“,i?“P] 


(A.l) 


where we have noted that AU comes from the metric per¬ 
turbation and are using the same notation as in Sec. IV 
of [1^ , where ~ denotes that we are neglecting any terms 
that do not lead to transcendentals and logarithms (in¬ 
cluding the overall scaling). We have suppressed the de¬ 
pendence of R™ and R’^p on i and m here (and in similar 
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expressions later), for simplicity. Note that in the ex¬ 
pressions in previous sections we denote i?'” and by 
Rh and i?oo, respectively. Also, we have (Eq. (166) in 
Sasaki and Tagoshi M) 

= K^R’^ + (A.2) 

{Ki, and Rq are given in, e.^Eqs. (6) and (7) in [s^) 
and [Eqs. (4.1) and (4.9) in [s^, evaluated for |s| = 2] 


i?"P = 


S^R’^ - 


where we have defined 


sin7r(i/ -|- ie) 
sin7r(i/ — ie) ’ 


(A.3) 


(A.4) 


discussion, where we are concerned with the simpli¬ 
fication. 

Now (recalling that e = 2mv^ and wro = mv), we have 


~ (1 - 2^2)e-—(2mz;)'^^%±i^. 

^ ^ ^ r(i-f2i/) 

(A.9) 

Thus, the RqRq'^ ^ term in Tim contributes 

^^,r(l + i.-ie)r(l-i. + ie) 


r(i-f 2^)r(i - 2z/) 


-T, 


(A.IO) 


where 


T := -4--+ c.c. . 


021771/ 


(A.ll) 


and (Eq. (23) in Sasaki and Tagoshi 0) 

R'^P] - (A.5) 

denotes the Wronskian of R™ and R'^^. Here [Egs. (157), 
(158), (168), and (170) in Sasaki and Tagoshi [Sg, noting 
that K = 1 for Schwarzschild] 


This cannot contribute any eulerlog terms (the expansion 
of the gamma functions does not contain ay), so we leave 
it alone. 

The (i?c)^ term in T^rn, on the other hand, does give 
exactly the eulerlog contribution found in T^. Specifi¬ 
cally, it gives 


^trans ^ ^ie 

(A.6a) 


B”® ~ {K, - 

(A.6b) 


7(1 + 27J 

■ ^ -^el2A^ut2 r(l + + U) 

^ ® r(i + i.-k)’ 

(A.6c) 

= X exp 

2i^ eulerlog^ (u) -I- 2'Kmv^ -I- ^ 

n—2 

A''_ - 

(A.6d) 




, 

n 

(A.12) 


We thus have 


where 


W[i?“, R'^P] - {K^ - A-'+A'i 


- {K, - ie-'-"5,if_,_i) 


r(l + l/ + le) - 7 re 
T{l + u- ie) 

(A.7) 


so we can write Eq. (lA.ll) as 

,,r(l + i. - ie) {S^R-^ - ie-‘^i?c"-') 




r(l-f ie) 

Rr''~ 

if. 


1 + 


1 — ie” 


^5, 


if-.-i 

if. 


Si, -l- 
-I- c.c. . 


o2i7ri/ 


(A.8) 


The K-^-iRq'' ^/{K^Rq) term is likely the ori¬ 
gin of the contribution to the 

simplification (just as it is for the V(,m simpli¬ 
fication in d^), since K-i,-iRq''~^ / {K^R'Q ~ 

(2u^)^)^ {gamma function terms} (cf. Eqs. (27a) and (27c) 
in [ 33 ). As these K-n-ij terms only contribute at 
higher orders, as discussed in Sec. IV of [s^, we shall 
thus omit the final fraction in the product in the ensuing 


- (l-2u2)-"™"%-2i™ 

v := _ _-f C.C., 


1 




,,3\n 


Q := {—V — 2imv'^)^ + {—u + 2\mv'^) 


2 (- 27 ”. 

(A.13) 


[Here we have abused notation in the “physicist’s way” 
with eulerlog,^ again, writing eulerlog^ (u) := 7 -|- 

log( 2 mt;), which is not what one would obtain when sub¬ 
stituting V for the argument of either of the previous two 
definitions, but, of course, agrees with them when one 
substitutes the values of i? or u corresponding to this 
u.] Unfortunately, the process of obtaining the full sim¬ 
plification from a study of the pieces entering the MST 
computation is obviously more subtle in this case than 
it is for the energy flux (discussed in [s^): The remain¬ 
ing terms in the expansion of this quantity [i.e., leaving 
off the factor] are not those found from a 

study of the expansion of and given in Eq. (fTHl) . The 
terms obtained from this expansion are more numerous 
and do not have the correct coefficients. The leading 
term indeed has the factor of 1 /j^, but none of the other 
terms seem to match. 
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